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We study ground-state properties of the two-dimensional Hubbard model at half filling by 
improving variational Monte Carlo method and by implementing quantum-number projection 
and multi-variable optimization. The improved variational wave function enables a highly ac- 
curate description of the Mott transition and strong fluctuations in metals. We clarify how 
anomalous metals appear near the first-order Mott transition. The double occupancy stays 
nearly constant as a function of the on-site Coulomb interaction in the metallic phase near the 
Mott transition in agreement with the previous unbiased results. This unconventional metal at 
half filling is stabilized by a formation of "electron-like pockets" coexisting with an arc struc- 
ture, which leads to a prominent differentiation of electrons in momentum space. An abrupt 
collapse of the "pocket" and "arc" drives the first-order Mott transition. 

KEYWORDS: variational Monte Carlo method, strongly correlated electron systems, Hubbard model, Mott 
transition, Fermi arc 



When electrons become localized by the strong 
Coulomb interaction, the Mott insulator appears after a 
metal-insulator transition called the Mott transition. 1 ' 2 
The Mott transition is found in many materials and sys- 
tems, such as transition metal oxides, 2 organic conduc- 
tors, 3 and 3 He layered systems. 4,5 

Filling-control Mott transitions at zero temperature 
have been studied by the auxiliary-field quantum Monte 
Carlo (AFQMC) method in the Hubbard model on a 
square lattice. 6, 7 The transition shows a continuous char- 
acter with divergences of the compressibility and the 
antiferromagnetic correlation length. Bandwidth-control 
Mott transitions have been studied by the path-integral 
renormalization group (PIRG) method. 8 ~ 10 The results 
have shown the first-order as well as continuous Mott 
transitions. The Mott transition has also been studied by 
the dynamical mean-field theory (DMFT). 11 The DMFT 
has shown that the Mott transition at the critical end 
point at nonzero temperature is consistent with the Ising 
universality class. 12 Extensions of the DMFT to include 
spatial correlations, 13 ~ 15 such as the cellular dynamical 
mean- field theory (CDMFT), have also suggested the ex- 
istence of both types of the first-order and continuous 
transitions. The Mott transition at finite temperatures 
has been studied by a phenomenological effective the- 
ory 16 and a mean-field theory as well. 17 ' 18 They have 
found unconventional criticalities at a marginal quantum 
critical point (MQCP), which arises from an interplay of 
topological transition and symmetry breaking. Experi- 
mentally measured critical exponents of an organic con- 
ductor, k-(BEDT-TTF) 2 Cu[N(CN) 2 ]C1 are in agreement 
with those of the MQCP. 19 

While crucial properties of the Mott transition, such 
as criticalities, have been elucidated by recent extensive 
studies, underlying electronic structure of metals near 
the Mott insulators has not fully been understood. This 
is an important issue to be solved, because competing or- 
ders and quantum phases near the Mott transition must 



be consequences of the underlying unconventional elec- 
tronic structure. Filling-control Mott transitions gener- 
ating the high-temperature superconductivity in the cop- 
per oxides, and bandwidth-control Mott transitions gen- 
erating the superconductivity and the quantum spin liq- 
uid in the organic conductors 20 are such consequences. 

Recently, electron differentiations in momentum space 
around the Mott transition are extensively studied both 
in experimental 21 and theoretical 22 ' 23 studies as a key 
to understand the unconventional electronic structure. 
The angle resolved photoemission spectroscopy has re- 
vealed a truncation of the Fermi surface called Fermi 
arc, in the underdoped region of the high-temperature 
copper oxide superconductors (HTSC). 21 This trunca- 
tion is a typical differentiation where electrons at differ- 
ent points of the Fermi surface start showing distinct be- 
havior. The origin of the Fermi arc and the pseudogap in 
the HTSC has been extensively studied by using the ex- 
tension of DMFT, such as the CDMFT. 24 - 25 In order to 
capture the differentiation, an accurate treatment of spa- 
tial correlations and high resolution in momentum space 
are crucially important, while extensions of CDMFT to 
larger spatial sizes are difficult. For this purpose, the vari- 
ational Monte Carlo (VMC) method 26 offers an alter- 
native advantageous approach, where large system sizes 
are tractable even with strong interactions and geomet- 
rical frustration effects. It offers high-resolution results 
in momentum space. However, the bias inherently and 
inevitably contained in the assumed variational form of 
the wave functions is a severe limitation in the VMC 
method. Reducing the biases is crucially important in 
studies with variational wave functions. 

In this paper, we study the electron differentiation of 
the Mott transition by improving the VMC method. We 
apply our improved variational wave function 27 to the 
two-dimensional Hubbard model with geometrical frus- 
tration effects. We find that highly fluctuating metals 
with large amplitude of double occupancy of electrons is 
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retained near the Mott transition. This unconventional 
metal at half filling is stabilized by the electron differ- 
entiation through a formation of "electron-like pockets" 
together with an arc structure. 

The Hubbard model on a square lattice is defined as 



0.06 



n = ~^2 ^Aa^a + Uj2 n iT 



nti, 



(1) 



where c i(T (cj ct ) is the creation (annihilation) operator on 



the i-th site with spin a and rii a 



W-T-Ci, 



is the number 



operator. The transfer integral is taken as iy = t for the 
nearest-neighbor sites and ty = t' for the next nearest- 
neighbor sites. We take N s = Lx L sites with the bound- 
ary condition periodic in the x direction and antiperiodic 
in the y direction (periodic-antiperiodic boundary condi- 
tion). Throughout this paper, we consider the half-filled 
case n = (1/N S ) £^(n l(J ) = 1 at t'/t = -0.3. 

The variational wave function used in this study has 
the following form: 27 

l^)=^d^G£ S=O |0pair}, (2) 

where Vq, ^d-h' an< ^ ^J are ^ ne Gutzwiller factor, 28 the 
doublon-holon correlation factor, 29,30 and the Jastrow 
factor, 31 respectively. These factors are defined as 
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where g, oti^., and wy are variational parameters. We as- 
sume that v^ = v(ri — rj) depends on the displacement 
r i — rj . Here, £\L\ is a many-body operator which is di- 
agonal in the real space representations. When a doublon 
(holon) exists at the i-th site and m holons (doublons) 
surround at the £-th nearest neighbor, C, \ returns 1. In 



;W 



i(m) 



other cases, £., 1 returns 0. The spin quantum-number 

projection C s=0 restores the SU(2) spin-rotational sym- 
metry and generates a state with total spin S = 0. 32 ' 33 
The one-body part |<^ pa i r ) is the generalized pairing wave 
function defined as 



E^ (1) ( fe )4Ai + E^ (2) ( fc ) 

.feeBZ fceAFBZ 



,t J 
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with the conditions tp^(—k) — tp^(jk) and <// 2 )(— fe) = 
tp( 2 \k). Here, BZ and AFBZ denote the Brillouin zone 
and the folded antiferromagnetic (AF) Brillouin zone, re- 
spectively. The pair orbitals (p^ (fe) and tp^ (fe) are vari- 
ational parameters. Since the PIRG method shows that 
the ordering vector of the AF insulating phase in the in- 
terval, < -t'/t < 0.5, is Q = (tt, tt), 8 - 10 the vector Q in 
eq. (6) is chosen as Q — (tt, tt) in this study. All the varia- 
tional parameters (g, a,K, Uy, (p^(k), and (p( 2 \k)) are 
simultaneously optimized by using the stochastic recon- 
figuration method. 34 The variational wave function \ip) in 
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Fig. 1. (Color online) (a) Squared staggered-magnetization m^, 
squared local moment w»? oca i, (b) double occupancy D, and 
charge gap A c in the thermodynamic limit as functions of the 
on-site interaction U/t at t'/t = —0.3. Each solid curves are fit- 
ted by the third-order Bezier curve. The inset is an enlargement 
of (b) for D near the transition. 

eq. (2) allows to describe various states such as paramag- 
netic metals (PM), antiferromagnetically ordered states, 
and superconducting states with any wavenumber (spa- 
tial) dependence of gap functions within a single func- 
tional form. Moreover, \ip) enables efficient treatment of 
quantum fluctuations with long-ranged as well as short- 
ranged correlations. For example, the peak value of the 
spin structure factor in the doped Hubbard model shows 
quantitative agreement with the results obtained from 
the unbiased AFQMC method. Detailed discussions and 
extensive benchmarks are reported elsewhere. 27 

In order to capture the metal-insulator and mag- 
netic phase transitions, we calculate the double occu- 
pancy D(N S ) = 0-/N s )^2Aniinii), spin structure fac- 



tor S(q,N s ) = (l/3N s )J2 g {Si-Sj 



%q-[Ti 



> , and charge 



? a P A C (N S ) = [u((N s + 1)/JV„) - M ((7V 8 - l)/JV„)]/2 
for Ar s -site systems, where Si is the spin operator and 
the chemical potential /i is given as /i((2M — 1)/N S ) — 
[E(M,M)-E(M-l,M-l)}/2. Here, E(N h Ni) is the 
variational energy with the number of up (down) spin 
JV-f (Ni). The double occupancy in the thermodynamic 
limit D is estimated by fitting the finite-size data in the 
form D — D(N S ) oc L~ l , because the PIRG results for the 
frustrated Hubbard model have succeeded in the extrap- 
olation by this form and estimating the critical point for 
the MIT. 8,9 The staggered magnetization in the thermo- 
dynamic limit m s is also estimated by fitting the scaling 
form m 2 ~S(Q, N s )/N s ex i" 1 with Q = (vr, tt) by follow- 
ing the scaling of the spin wave theory if the long-range 
order is present. 35 We also calculate the local moment 
Wlocal defined as mi OC aJ = \/JSfSf) = (1/2)V«- 2L>. 
For the charge gap in the thermodynamic limit A c , we 
use the scaling function A c — A C (N S ) ex L^ 1 as in the 
AF Hartree-Fock gap equation. 6 The extrapolation to 
the thermodynamic limit is performed by using several 
choices of system sizes up to L = 16. 

Figure 1 shows the U/t dependence of D, m 2 , m 2 ocal , 
and A c . The first-order metal- insulator transition takes 
place at U c /t = 3.3 ±0.1. The magnetic phase transition 
takes place at the same critical point within our reso- 
lution. The nonmagnetic insulating phase is not clearly 
identified in our variational results. Although mi oca i 
gradually increases as a function of U/t, the growth is 
strongly suppressed in the metallic phase near the metal- 
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insulator transition. In the metallic region, S(Q, N a ) at 
Q = (vr, 7r) rapidly grows when U/t approaches U c /t, al- 
though the true long-ranged order is not achieved within 
metals. This growth is, however, to a large extent, as- 
cribed to the growth in the longer ranged part in the 
real space correlation, while the shorter ranged part of 
antiferromagnetic correlations does not grow equally. In 
this overall tendency, the local moment ™f oca i, that is 
the shortest-ranged part, shows nearly flat dependence 
on U/t. This seems to optimize and reconcile in gain- 
ing the kinetic energy by keeping the electron-pocket- 
like and arc-like structure of the Fermi surface as we 
describe below. In other words, the flat U/t dependence 
of m Lai = hmAr s ^oo(l/AyE,j 5 '(<Z; A 's) with growing 
S(Q, N s ) means a compensating reduction of S(q, N s ) 
at q 7^ Q. This reduction suppresses the electron scat- 
tering by the spin fluctuations and keeps the coherence 
of the carrier, leading to the gain in the kinetic energy. 

Figure 1 shows remarkable quantitative agreement 
with the unbiased results obtained from the PIRG 
method. 8 Although the value t'/t in the PIRG data in ref. 
8 is not the same as that of the present results, the values 
at t'/t = —0.3 are rather precisely interpolated from the 
results at t'/t = -0.2 and t'/t = -0.5 given by the PIRG 
method. For example, U c /t estimated by interpolation of 
the PIRG results suggests U c /t ~ 3.6 ± 0.1, which may 
be compared with the present estimate U c /t = 3.3 ± 0.1. 
The absence or only tiny interval of the nonmagnetic 
insulating region at t'/t = —0.3 is also consistent each 
other. Furthermore, D near U c /t in the metallic phase 
stays nearly constant around D ~ 0.2 over an interval 
2 < U/t < 3.3 contrary to a naive expectation. This flat 
U/t dependence of D also remarkably agrees with the 
PIRG results. In addition, Fig. 1 shows convex growth of 
the squared staggered-magnetization mj? from U c /t for 
U/t > U c /t. These features have been observed by the 
PIRG method at t'/t = —0.2 as well. In our results, 
the squared staggered-magnetization at U/t = 4.0 and 
t'/t = -0.3 is ml = 0.020 ± 0.002. In the PIRG results, 
the value at U/t = 4.0 and t'/t = -0.2 is ml ~ 0.025. 
When we consider the difference of t'/t, this is a remark- 
able quantitative agreement. Our results suggest that our 
variational wave function enables quantitatively accurate 
descriptions in the ground state. 

On the other hand, in the previous VMC calculation 
for this model, the critical value of the Mott transition is 
much larger (U c /t ~ 6.7) . 36 In addition, the double occu- 
pancy in ref. 36 linearly decays to the transition point as 
a function of U/t with a substantial slope (dD/d(U/t) ~ 
—0.033). The variational wave functions employed in the 
literature include many-body correlations only in much 
restricted forms, such as the short-ranged doublon-holon 
factor. Such a restricted form does not sufficiently take 
into account fluctuations, which are strongly enhanced 
around the Mott transition. By introducing a large num- 
ber of variational parameters that scales linearly with 
the system size, the Gutzwiller-Jastrow factor as well as 
the one-body part allow much more accurate treatment 
of fluctuations . DMFT results also show large U c /t ~ 11 
even for t'/t = and fail in capturing the plateau of D 
at U/t < U c /t n Although the CDMFT takes into ac- 




Fig. 2. (Color online) (a) U/t dependence of momentum distri- 
bution n(fc) for L = 14 and t'/t = —0.3. Error bars are com- 
parable to the symbol sizes. The filled points are values of 
n(k) in eq. (7) and n(k) obtained by the linear interpolation 
of two nearest fc-points along the vertical direction of each _T- 
M-X-T line. PM and AFI denote the paramagnetic metal and 
the antiferromagnetic insulator, respectively. Dashed curve in 
the inset is the Fermi surface for t'/t = —0.3 and U/t = 0. 
The three-dimensional plot of n.(fc) at (b) U/t = 3.2 and (c) 
U/t = 3.4 for L = 14 and t'/t = -0.3. The value n(k) at 
{k = ({2£ x - l)ir/L,(2£ y - l)vr/L), {2£ x n/L, 2£ y n/L) \ £ x ,£ y = 
—L + 1, — L + 2, ■ ■ ■ ,L — 1, L} is obtained from the linear inter- 
polation of n(fc) in eq. (7). 

count spatial correlations to some extent, the plateau of 
D around U c /t has not been captured yet in the CDMFT 
studies. 37,38 The exact diagonalization study 39 does not 
capture this behavior either because of the limitation of 
the system size. Sufficiently large system size over the 
correlation length of the fluctuations is important for re- 
producing the sufficient coherence and the plateau of D. 
In order to analyze the electron differentiation in 
this region, we calculate the momentum distribution 
n{k) = {l/2N s )J2 hJ Jclc Jr7 }e lk < r ^ r ^. We assume 
that n(k) has fourfold symmetric structure. Since we 
calculate n(k) under the periodic-antiperiodic boundary 
condition, this assumption allows twice as many fc-points 
in the BZ for the L x L system: 

M ■(*■§■ 

L + 1, -L + 2, • • • , L - 1, L. The VMC 

results for L = 14 are shown in Fig. 2. The momentum 
distribution n(fc) shows a characteristic behavior near 
the Mott transition. In the metallic phase, the amplitude 
of n(fc) around M points (fc = (±7r, 0), (0, ±7r)) is kept 
large, generating an "electron pocket-like" structure just 
before the Mott transition as is seen in Figs. 2(a) and (b). 
Figures 3(a)-(d) are contour plots of |Vn(fc)|. We can 
estimate the Fermi surface and the coherence of elec- 
trons on the Fermi surface by the amplitude of |Vn(fc)|. 
The contour plots remarkably show the "arc-like" struc- 
ture in the metallic phase close to the Mott transi- 
tion (Fig. 3(c)). Although our system size is not suf- 
ficient, the position of the Fermi surface looks rigid 
and very similar to that of the non-interacting sys- 
tem, contrary to the "deformation to the nesting" 
obtained in the renormalization-group method in the 
weak-coupling regime. 40 The rigidity enhances the co- 
existence of "pocket" and "arc." The "arc" around 
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(a) U/t = 1.0 
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(b) U/t = 3.0 
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Fig. 3. (Color online) Contour plots of |Vn(fe)| for L = 14 and 
£'/£ = —0.3. Here, | Vn(fe)| is calculated by using the value from 
the bicubic interpolation of n(fc) in eq. (7). Dashed curves are 
the Fermi surface for t' /t = —0.3 and U/t = 0. 

(±7r/2, ±7r/2) coexisting with the "electron pocket-like" 
structure around (±7r, 0) and (0, ±ir) discussed before 
is retained near the Mott transition. The coexistence 
implies that, although the long-range magnetic order 
is absent, the gap structure of the charge excitation 
already gets formed even in the metallic phase as a 
precursor of the Mott gap. As in semimetals, the pre- 
formed Mott gap generates electron and hole pockets. 18 
The arc structure emerges as a part of the hole pock- 
ets around (±n/2, ±7r/2) where outer half of the pocket 
is lost presumably because of a strong damping. The 
coexistence of the "electron pocket" and "arc" directly 
causes largely retained D and suppressed ?7ii oca i in Fig. 
1. This "semimctallic pocket" structure creates electrons 
and holes in the upper and lower Hubbard bands, respec- 
tively. As a result, excitons (or doublon-holon pairs) are 
generated in the background of the Mott insulator, where 
D is largely retained. The first-order Mott transition ap- 
pears as a sudden collapse of the "semimetallic pocket" 
structure. The similar abrupt change of n(k) is also 
seen in the PIRG results. 8 The Mott transition emerging 
as shrinkage of electron-hole (or doublon-holon) pockets 
supports the topological character of the transition pro- 
posed in ref. 18. 

In summary, we have studied the electron differenti- 
ation in momentum space around the Mott transition 
by using the improved VMC method. Our variational 
results show quantitative agreement with the unbiased 
method. Especially, our variational wave function enables 
to treat short-ranged as well as long-ranged spin/charge 
fluctuations, which are strongly enhanced near the Mott 
transition. As a result, we have captured the flat U/t 
dependence of double occupancy D and abrupt change 
of momentum distribution n(k) at the first-order Mott 
transition. The momentum distribution shows "electron 
pocket-like" structure around k = (±7r, 0), (0, ±7r). The 
arc structure appears around k = (±7r/2,±7r/2) in 
|Vn(fe)| in the metallic phase near the Mott transition. 
The abrupt collapse of these parts drives the first-order 
Mott transition. The coexistence of "electron pocket- 
like" structure and "hole-like arcs" with the retained 
plateau of the double occupancy D is the key aspects 
of the electron differentiation in momentum space. Clar- 
ifying the relation between the above differentiation and 
instability toward superconductivity is one of the most 



important issues left for future studies. 
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